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Abstract. The dynamics of coUisionless galaxy can be described by the Vlasov- 
Poisson system. By the Jean's theorem, all the spherically symmetric steady 
galaxy models are given by a distribution of <&(-©, L), where E is the particle 
energy and L the angular momentum. In a celebrated Doremus-Feix-Baumann 
Theorem [7], the galaxy model <^{E,L) is stable if the distribution 3> is mono- 
tonically decreasing with respect to the particle energy E. On the other hand, 
the stability of ^{E, L) remains largely open otherwise. Based on a recent ab- 
stract instability criterion of Guo-Lin we constuct examples of unstable 
galaxy models of f{E, L) and / {E) in which / fails to be monotone in E. 



1. Introduction 

A galaxy is an ensemble of billions of stars, which interact by the gravitational 
field they create collectively. For galaxies, the collisional relaxation time is much 
longer than the age of the universe ([5]). The collisions can therefore be ignored and 
the galactic dynamics is well described by the Vlasov - Poisson system (coUisionless 
Boltzmann equation) 



(1) dtf 
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where {x,v) e M"^ x M'^, f{t,x,v) is the distribution function and Uf{t,x) is its 
self-consistent gravitational potential. The Vlasov-Poisson system can also be used 
to describe the dynamics of globular clusters over their period of orbital revolutions 
([9]). One of the central questions in such galactic problems, which has attracted 
considerable attention in the astrophysics literature, of [5], [6], [9], [17] and the 
references there, is to determine dynamical stability of steady galaxy models. Sta- 
bility study can be used to test a proposed configuration as a model for a real stellar 
system. On the other hand, instabilities of steady galaxy models can be used to 
explain some of the striking irregularities of galaxies, such as spiral arms as arising 
from the instability of an initially featureless galaxy disk ([5]), Ql8j). 

In this article, we consider stability of spherical galaxies, which are the simplest 
elliptical galaxy models. Though most elliptical galaxies are known to be non- 
spherical, the study of instability and dynamical evolution of spherical galaxies 
could be useful to understand more complicated and practical galaxy models. By 
Jeans's Theorem, a steady spherical galaxy is of the form f{x, v) = fi{E, L), where 
the particle energy and total momentum are E = \\v\'^ + U{x), L = \x x v\ , and 
Uf_i{x) — U {\x\) satisfies the self-consistent nonlinear Poisson equation 



Af7 = 47r / fi{E,L)dv 
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The isotropic models take the form f{x,v) = IJ,{E). The case when ^'{E) < (on 
the support of n{E)) has been widely studied and these models are known to be 
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Stable ([T], [2], [2], [2], [8], [20], [14], [19], [13]). To understand such a stability, we 
expand the well-known Casimir-Energy functional (as a Liapunov functional) 

(2) ^ ^ / / ^^^^ + ^ / / '"''^ " ^ / 
which is conserved for all time with 

(3) Q'it^iE)) ^ ~E 

(this is possible only if fi' < 0!). Upon using Taylor expansion around a steady 
galaxy model of [f{E, L),U], the first order variation vanishes due to the choice of 
([3]), and the second order variation takes the form of 

The remarkable feature for stability lies in the fact that ^^[3,3] > if < 
(T], f2], [7], [2]). This crucial observation leads to the conclusion that galaxy 
models with monotonically decreasing energy /i' < are linearly stable. On the 
other hand, for the case fj, is not monotone in E, ^ breaks down, and H'^ is 
not well-defined, which indicates possible formation of instability. It is important 
to note that a negative direction of the quadratic form ^^'[5,5] < does not 
imply instability. In 112', an oscillatory instability was found for certain generalized 

polytropes / {E,L) = /oi"^" (-Bo - Ef^ with n < |,m < 0, by using N-body 
code. This instability was later reanalyzed by more sophisticated N-body code in 
[3]. Despite progresses made over the years (e.g., [3], [10], [17], [12), no explicit 
example of isotropic galaxy model fi{E) are known to be unstable. 

The difficulty of finding instability lies in the complexity of the linearized Vlasov- 
Poisson system around a spatially non-homogeneous fi{E, L) : 

(5) dtg + v-V^g-V^Ug-V^n~V:cU -V^g^Q, AUg = in g{x,v)dv 

for which the construction for dispersion relation for a growing mode is mathemat- 
ically very challenging. In a recent paper ([11]), a sufficient condition was derived 
rigorously as follows: 

Theorem 1. Let [f{E^L)^U] be a steady galaxy model. Assume that f{E,L) has 
a compact support in x and v, and /x' is bounded. Define auxiliary quadratic form 
for a spherically symmetric function ipilxl) "-^ 
(6) 

where ri{E,L) and r2{E,L) are two distinct roots to the equation 
(7) E - U{r) - L^l2r^ = 

and the average (j) is defined as 

rr2{E,L) (f>{r)dr 

- Jri{E,L) ^2(E~Uo{r)~Ly2r^) 

Jri(E,L) ^2(_E-(7o(r)-L2/2r2) 

// there exists 4>{\x\) such that [Ao<f>^(j)] < 0, then there exists an exponentially 
growing mode to the linearized Vlasov-Poisson system l^i- 
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The quadratic form [Aqc/), </>] in the above instabihty criterion involves dehcate 
integrals along particle paths. The purpose of the current paper is to use numerical 
computations to construct two explicit examples of galaxy models for which a test 
function (j) satisfying [Ao(/), 0] < exists. This ensures their radial instability by 
Guo-Lin's Theorem [TJ The first example is an anisotropic model with radial insta- 
bility. There are two differences to the oscillatory instability found in [T? and [3]: 
here the distribution function is non-singular and the instability is non-oscillatory. 
The second example is a non-singular isotropic model with radial instability. To 
our knowledge, this provides the first example of unstable isotropic model. 

Even though the galaxy models studied here are not actual ones observed, our re- 
sults demonstrates how to apply Guo-Lin's Theorem [1] to detect possible instability 
for a given galaxy model. It is our hope to foster interactions between mathemat- 
ical and astronomical communities and to advance the study of instability of real 
galaxy models. 



2. Examples of Unstable galaxy models 

Example 1. Unstable Galaxy Model Depending on E and L. We define 
the distribution function /o {E, L) ~ fi{E)L'^ where 

{0 E <4: 

2.25(i?-4)2 4<E<4A 
{b-Ef 4.4<i?<5 
S > 5 

The graph of ^m{E) is showed in Figured] below. Choosing U{Q) — 3, we numerically 



0.4 




Figure 1. ^j.[E) 



calculate U{r), and the graph of U{r) is showed in Figure [2] By choosing the test 
function — e^^, then numerical computation below shows 

[Ao0, 0] < TT - 45 < 0, 

and hence /o is unstable by Theorem [T] 



4 



ZHIYU WANG, YAN GUO, ZHIWU LIN, AND PINGWEN ZHANG 




Figure 2. U{r) 

Example 2. Unstable Galaxy Model Depending Only on E. Let fo{E) = 
^{E), where 



(9) ^iiE) 



E <a 

Cit\{E~af a<E<Ei 
Ci{Eo-Ef Ei<E<Eo 
E> Eo 



We choose 

Ci = 2, Eo = 5.1, a = 1.9, k = 2.01, h ^ 1.5. 
Choose Uq = —15.1, the graphs of n{E) and U are shown in Figures |3] and HI We 




Figure 3. fi{E) 

choose the test function 
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Figure 4. U{r) 



with wi, W2, ...w? given by (fT6|). Then [Ao(/ 
profile follows from Theorem [TJ 



-24. 016 < 0, and instability of the 



3. Numerical implementation 

3.1. Numerical computation for Example 1. By the resuhs in [21], for a steady 
state with /o {E, L) ^ fi {E) X^' of Vlasov-Poisson system, the steady potential 
U (r) satisfies the equation 



(10) 



U'{r) = 



9i+i/2{U{s))ds, r>0 



with gm(u) = fi{E){E — u)™-dE, u G (— oo, oo), and 



Ca,fc 



Jo 



r(a + i)r(b + i) 

r(a + 6 + 2) ' 



a > -1, b > -1, 



where F denotes the gamma function. The boundary condition is limj._).oo U (r) = 0. 

The computation of finding (/> with [Aq4>, 0] < is carried out in the following 
steps. 

Step 1. Computation of potential U. For Example 1, /o (-E, L) = ^ (E) where 
fi (E) is given by ([8]) , so Z = 2 and 



m+i/2{u) 



4.4 1-5 



2.25{E-4:y{E-uy-''dE 



{b~ EYiE-uY-HE. 



Letting M = 2'+'^/27r2Q._i/2, we obtain from ([TOl) 

r''U"{r)+2rU{r) = Mr2'+25;+i/2(;7(r)), 
which is equivalent to 

r ' 
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Table 1. Lmax for different E 



E 


Lfnax 


3.2 


0.452787235493472 


3.4 


0.731586885359610 


3.6 


0.969996239842229 


3.8 


1.186531828839693 


4.0 


1.388987957485868 


4.2 


1.581703995145215 


4.4 


1.767449993848502 


4.6 


1.948250754546466 


4.8 


2.125722162034805 


5.0 


2.301341325775183 



We use the boundary conditions U'{Q) — and U{Q) — 3. Let ui{r) — U{r), U2{r) 
u'l (r) , we transform the above equation to the following system 



(11) 



Kir) 
Kir) 



U2 



U2 + Mr^^gi+i/2{ui) 



with Mi(r) = 3 and M2(0) = 0. We remark that we can subtract the finite limit of 
U at infinity from U and redefine Eq accordingly. We apply Runge-Kutta method 
to solve equations PT|) . as follows 



h 
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Vn+l 

K2 = fiXn + 

Ki = f{xn + h, y 



{Kl + 2K2 + 2K3 + Ki) 



2 ' 
2'' 



-1^1) 
--2K2) 
hK2). 



First let h = 0.1 to get the values of U and U' at points Xn ~ nh, then we use 
piecewise cubic Hermite interpolation to get an approximation of U (r). 
Step 2. Computation of roots of the equation 



(12) 



E-U{r) -L^jlr^ = 0. 



For fixed E^ this equation has two solutions r\ < r2 when L < imax, one solution 
r* when L = Lmax and no solution when L > L„iax. Here, Lmax iE) — y/r*^U'{r*) 
and r* (E) satisfies the equation 



(13) 



E-U{r) - -U'{r) = 0, 



which comes from the combination of the equations —U'{r) + 1? jr^ = and ([T2|) . 
We employ Newton method to find the unique root r* of by the Newton 
iteration 

E - U{r) - W'{r) 



'„+! 



\U'{r) - W"{r) 



Chose an initial point 7'o = 1, we get r*(with the stopping criterion |r„ 
Table [Hand Figure [5] show the relation of Lmax to E. 



< 10- 
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Figure 5. relationship of Lmax and E 




Figure 6. E=A.5,L=0.5 Figure 7. E=i.8,L=1.2 

For L < L„iaxi we apply Newton's iteration 

E - U{r) - l?l2r'^ 



rn+i 



-U'{r) + L^/r^ 

to solve (|T2|) with the choice of initial value 

_ J 0.1 when L > 0.05, 

~ \ 0.001 when L < 0.05 ' ^^'^ ~ ' 

and the stopping criterion |r„ — r„+i | < 10"^". We show two graphs of E ~ U{r) — 
I? jlr^ in Figures [6l and [7] and the computation results are 

ri = 0.288701795314679, ra = 1.843722113067511 for Figurel 

ri = 0.634700793130700, = 1.938943620330099 for Figure [7] 

Figures [8] and [9] show how n and {r\ < r2) change with respect to L, when E 
is given. (Note that the equation ()12p has a unique root r* when L = Lmax, thus 
when L approaches Lmax, the distance r2 — ri tends to zero.) 
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Figure 8. £'=3.4, Figure 9. E^A, 

how ri , r2 change with L how ri , r2 change with L 

Table 2. the value of I{E,L) for given E and L 



{E,L) 


inter grand{E, L) 


(4.1,0.6) 


0.002530474111772 


(4.3,1) 


0.053619816035575 


(4.5,1.2) 


-0.086444987713539 


(4.7,1.3) 


-0.086591041429108 


(4.9,1.7) 


-0.066267498867856 



Step 3. Computing A{(j),(j}) in Denote the integrand 
I{E,L)^UE,L) ' ' 



iri[EX) ^J2{E ~ C/o(r) - L2/2r-2) 

Choose (j) = e^*". Table [2] shows some of the values oi I{E, L). 

For given E, the integration interval oiI{E,L) in ([6]) is [0,Lma2;] for L. When 
L = and L = Lmax, has only one solution instead of two. To avoid this 
problem, we integrate I(£, L) for L in the truncated interval [ai , Lmax ~ ^2] , where 
ai,a2 > 0. Choose 02 = 10^'^, and compare the value of °" I{E, L)dL when 

fli — 0.01 and ai = 0.001 in Tabled We see that the error is approximately 10^*". 

Table 3. Value of J^"'"^^'"'' I{E,L)dL for ai = 0.01 and oi = 0.001 





ai = 0.01 


ai = 0.001 


inter grandE{4:.3) 


0.057831884696931 


0.057831864357174 


inter grandE{4:.5) 


-0.084810940209437 


-0.084810952811088 


inter grandE{A.7) 


-0.095869703408813 


-0.095869701098749 


inter grandE{4:.9) 


-0.058302133990229 


-0.058302146458761 



Step 4- Error Estimate. Next we give a theoretical estimate for the error intro- 



duced by the truncation. Let A{E) — 



r{E) 



dr 



y/2(E-U{r)) 



, where r{E) is the unique 
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solution oi E — U{r) — 0. For 0(r) = e the error term for L on [0, ai] is 

.r2{E,L) 2LdrdLdE 



(14) 



3271^ / / ^' (E) / 
Ji Jo Jn 



64 
< — 
- 6 



riiE,L) y^2{E - U{r) - /2r^) 

TT^ \fi'{E)\ max { A {E),T{E,ai)}dE af. 

In the above, we use the estimate 

rr2{E.L) , 

nE,L) = / < max{A{ElT{E,a,)} 

due to the facts that hmi^o+ T {E, L) — A (E) and T {E, L) is monotone for L on 
the small interval [0, Oi]. We numerically compute the coefficient on the right side 
of ([H]) by choosing ai = 0.01, then 

^TT^ / |//'(£;)| max {A (£;) ,r(i;,ai)}(i£; = 3.787856674901141 X 10^. 
6 Ji 

Thus ai = 0.01 is enough to make the first error to be of order 10~^°. 
The error on [Lmax - a2,Lmax] is 



(15) 



. , /-^""^ . 4 n'^'^'^K -.0 2LdrdLdE 



L„ 



ri{E,L) y^2{E-Uir)-Ly2r^) 



< GAtt-' / L„,aAEr\^^'{E)\m^ix{^====,TiE,L,,,,^-a2)}dE a2. 

Ji [ V0eft(''*; Anax) J 

where the effective potential 

0eff(?-) = t/o(r)+LV2?-' 

Since 



T{E,L) < max < ,T{E,L^a^^ - 02) > , L G [Lmax - 02, L 

I VCff(^*;^max) I 



maxj ; 



due to the facts that lim^^L _ T (E, L) = . ^ and T (E, L) is mono- 

tone for L on [L^ax - 02, Lmax]- 

Numerically computing the coefficient before 02 in the right hand side of (|15p. 
we get 0.605275913474830. Choose 02 = 10-^ the right side of ^ if of order 
10~^, which is good enough for what we need later on. 

Step 5. Conclusion. We choose (j){r) = e~^. The first term of © is / |V0p = 
TT, and the second term is computed to be —31.733535998660550, by choosing 
h = 0.1 in Runga-Kutta Method when solving U{r). Therefore {A<j>,(j)) = ir — 
31.733535998660550 < 0. We gradually decrease h, and repeat the whole process 
to calculate the second term of ([6]) . The result is shown in Table HI Thus when 
h becomes smaller, the second term of ([6|) is less than —45. Combining with the 
error estimates of (fT4|) and (fTSj) in Step 4, we can ensure {Ao(p, (f>) < 0. 
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Table 4. results for different h 



h 


the second term 


0.1 


-31.733535998660550 


0.05 


-38.227746546049751 


0.025 


-41.848643231095494 


0.01 


-44.146480909452201 


0.005 


-44.933896897054382 


0.0025 


-45.331673900274509 


0.001 


-45.571655158578380 



3.2. Numerical Computation for Example 2. With profile ([9]), we choose (f>{r) 
to be a linear combination of several functions e^"''", 1 < « < n, (/<(r) = ^ Oie""'*", 
with Gi to be determined. In the quadratic form (|6]), the first term is computed to 
be 

and the second term is 

r2{E,L) 2LdrdEdL 



3271^ / f'^{E, L) / (V - V a,</.,) 

J JrKE.L) 



ri(£,L) ^ ^ - ' ^2{E - Uo{r) - Ly2r^) 

r2{E,L) ^ _ 2LdrdEdL 2 



V327rM /^(ii;,L) / (0, 

J JrKE.L) 



rUE.L) ^2{E - Uo{r) - Ly2r^) ' 
„ /■ , rr2{E,L) 2LdrdEdL 
J Jri(E,L) y^2{E -Uo[r) - l2r^) 

= E Cifl? + E '^'^ijO-iaj 

Now (Ao(/), 0) = + + -|- Cij)aiaj. The corresponding matrix for 

this quadratic form is 5 = (sy), where sa = bi + Ci, Sij — bij -\- Cij. Denote Amin 
to be the minimum eigenvalue of the matrix S, then Amin < guarantees that there 
exists (p such that (AqcI), (f) < 0. We choose 

In ([9|), there are five free parameters in [i: a,k,ti, Eo,Ci and Ei is determined 
by El = [Eq + tia)/{ti -\- 1). Using Uq (the initial value of the potential) as an 
additional parameter, there are totally six parameters in our calculations. The idea 
is to view Xmin as a function of these parameters Amm = Ami„(a, fc, ii, i?o, Ci, Uq). 
Our goal now is to find proper parameters such that \min < 0. The numerical 
methods are the same as in Step 1-Step 3 in the computation of Example 1, except 
that now we use self-adaptive Runge-Kutta Method to solve (fTT|) . We next verify 
that the following choice of parameters 

Ci =2, Ea^ 5.1, a = 1.9, k = 2.01, ti = 1.5, Ua = -15.1, 

would lead to Xmin < and {Aocj), cj)) < 0. In Table [5] we show the results of Amin 
under different accuracies, where ai and 02 are the parameters used to truncate 
[0, Lmax] to [ai,Lmax — ^2] ■ Wc also Calculate the eigenvector corresponding to 
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Table 5. Xmin under different accuracies 



integration accuracy 


ai 


02 


accuracy of tlic roots 




io-« 


10"^ 


io-« 


10-1" 


-1.808068979998836 x 10-" 


io-« 


10-^ 


10-^ 


lo-iiJ 


-3.403231868973289 x 10-" 


io-« 


10-^ 


io-« 


10-" 


-3.739973399333609 x 10-" 




10-^ 


10-*' 


10-" 


-3.196606813186484 X 10-" 


10-" 


10-^ 


io-« 


10-" 


-3.025099283901562 x 10-" 


10-^^ 


10-" 


10-^ 


10-" 


-3.087487453314562 x 10-" 


10-^^ 


10-*^ 


10-^ 


10-" 


-3.080483943121587 X 10-" 


10-^^ 


10-*^ 


10-^ 


10-" 


-3.089704058599446 x 10-" 


10-14 


10-'^ 


10-' 


10-" 


-3.089187445892979 x 10-" 



7p 
6 - 




D 0.6 1 1 5 2 2 5 3 3.5 4 4.5 5 



Figure 10. 0(r) 

Amin hy v= {vi,V2,V3,V4,vz,ve,V7)'^ , wliere 

vi = 41.702767064740 

V2 = -14.949618378683 
= 4.856846351504 
(16) V4 = -201.361293458803 

V5 = 571.694416252419 

ve = -723.292461038838 

V7 = 327.842857021134. 

Let 

We draw a picture of this function </> in Figure [Tni Using ^ , we obtain 

(A(/),<?!)) = 35.231484866282720- 59.247837621248109 < 0, 

when the integration accuracy is IQ-^'^, oi = 10-", 02 = 10-^, and the root accu- 
racy is 10-". 
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Step 4- Error Estimate. We obtain theoretical error estimate as in Example 1. 
For any test function 0, the error for L on [0,ai] is 



327r3 / n'{E) 



ai .r^iE.L) 2LdrdLdE 



{E\t,'{E)>0} Jo Jrt(E,L) V^i^ - - L'^ /2r'^) 

<327r^|max(/) — min(/)pa^ / ^' [E)'n\ayi{A{E)^T[E,ai)}dE 

J{E\fi'(E)>0} 

Numerically compute the right hand side, and note 

I max./) - mine/) I < 2niax|0| < 2^\v,\ = 3.7714 x 10^. 

Choose oi = 10^^, the error is no more than 

4.322 X 10^ X (3.7714 x lO^)^ x (lO"^)^ = 6.1474 x 10"^ 

Next we estimate the error for L on [Lmax — 02 , -^max] ■ When re [ri , r2] , we use 
|</) — 01 < max|0'|(r2 — ri) by the Mean Value Theorem. So the error is 

. f /■^™"- fr2(E^L) 2LdrdLdE 



{E\f^'{E)>0} JL^^,~a2 Jri{E,L) V^i^ - Uq - L^/2r^) 

< 64^3^2 / (E) i,„ax (E) {r2{E, L,„ax - 02) - ri{E, L„,ax - a2)f ■ 

J{E\ ^i'(E)>0} 

\ / ,,, ,^ r =^''^ -^max - 02) > dE . 
[ V/Cff(^*;imax) J 

Note that 

max < 5^ Iw,! = 9.4285 x 10^. 
Numerically computation of this error with 02 — lO"'' yields a bound of 

1.41722 X IQ-^^ X (9.4285 x 10^)^ = 1.3362 x 10"*^ 
So the total error is at most of order 10^^, which guarantees (Acj), (f)) < 0. 

4. Summary 

We study the instability of spherical galaxy models in the Vlasov theory for 
coUisionless stars. Based on the instability criterion of Theorem [1] and careful 
numerical computations, we have constructed two explicit unstable galaxy models 
fo{E,L) and fo{E) with distributions ^ and © respectively. In particular, fo{E) 
in Example 2 provides the first example of unstable isotropic galaxy which has 
not been found in literature. The instability in these examples are radial and 
non-oscillatory. Compared with the usual N-body codes of finding instability of 
galaxy models, our method only requires numerical evaluation of certain explicit 
integrals. Therefore, it is much more reliable and easier to implement. It is hoped 
that Theorem [1] can be employed to detest instability for other galaxy models in 
the future. 
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